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ABSTRACT 

Solid particles in protoplanetary disks that are sufficiently super-solar in 
metallicity overcome turbulence generated by vertical shear to gravitationally 
condense into planetesimals. Super-solar metallicities result if solid particles pile 
up as they migrate starward as a result of aerodynamic drag. Previous analyses 
of aerodynamic drift rates that account for mean flow differences between gas 
and particles yield particle pile-ups. We improve on these studies not only by 
accounting for the collective inertia of solids relative to that of gas, but also by 
including the transport of angular momentum by turbulent stresses within the 
particle layer. These turbulent stresses are derived in a physically self-consistent 
manner from the structure of marginally Kelvin- Helmholtz turbulent flows. They 
are not calculated using the usual plate drag formulae, whose use we explain is 
inappropriate. Accounting for the relative inertia of solids to gas retards, but 
does not prevent, particle pile-ups, and generates more spatially extended re- 
gions of metal enrichment. Turbulent transport hastens pile-ups. We conclude 
that particle pile-up is a robust outcome in sufficiently passive protoplanetary 
disks. Connections to observations of circumstellar disks, including the Kuiper 
Belt, and the architectures of planetary systems are made. 

Subject headings: planetary systems: protoplanetary disks — planets and satel- 
lites: formation — turbulence — methods: numerical 



1. Introduction 

The formation of planetesimals by gravitational instability within dense particle layers 
in circumstellar disks is an appealing alternative to collisional sticking of dust aggregates 
(Youdin & Shu 2002, hereafter YS02, and references therein). In one fell swoop, gravity 
promises to assemble sand-sized or smaller solids into kilometer-sized agglomerates. 

For gravitational instability to operate, two main requirements must be met. First, the 
disk must be sufficiently passive; turbulent velocity fluctuations in disk gas must be small 
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enough to allow solids to gravitationally settle toward the midplane. Kelvin-Helmholtz tur- 
bulence generated by vertical shear within the stratified particle layer threatens to violate 
this requirement (Cuzzi et al. 1993; Weidenschilling 1995). Such turbulence, however, can 
entrain only a finite amount of solids. Thus, a second requirement for gravitational insta- 
bility, first explored by Sekiya (1998), is that the surface density of solids relative to that 
of gas, Sp/Sg, must exceed a critical threshold, £ PiCr /£ g . Within the minimum-mass solar 
nebula (MMSN), £ PiCr /£ g lies ~5 — 20 times above the solar value of X p /X g ~ 5 x 10~ 3 . 

Many possibilities exist for achieving this enhancement of metallicity, including the 
return of condensed material to an accretion disk from a bipolar outflow and removal of 
dust-depleted gas in disk surface layers by magnetic accretion, photoevaporation, or en- 
trapment within a stellar wind. YS02 discuss the origins and merits of these possibilities. 
They highlight a promising and natural enrichment mechanism that arises from simple aero- 
dynamic drag. Solid particles rotate about the central star faster than does ambient gas 
because the latter is more sensitive to pressure gradients that are (usually) directed radially 
outward. This difference in mean flow velocities causes particles to be frictionally dragged 
by gas. Such drag induces orbital decay. By applying the Epstein aerodynamic drag law to 
millimeter-sized solids in the MMSN, YS02 conclude that as particles drift relative to gas to- 
ward the central star, they pile up: local enhancements of more than an order-of-magnitude 
in the particle surface density (metallicity) occur at small stellocentric radii. 

The Epstein drag law applies when relative velocities, i> re i, between a solid particle and 
surrounding gas are less than the gas sound speed, c g , and when the particle radius, a, 
is smaller than the gaseous mean free path, A. These conditions are satisfied for particles 
having a < 1 (r/AU) 2 75 cm in the MMSN, where r is the stellocentric distance. The Epstein 
drag force on a spherical particle reads 



where p g is the mass density of gas. YS02 find that millimeter-sized particles pile up to the 
point where their density exceeds the Roche density in ~10 5 yr, fast enough to occur within 
disk lifetimes (~10 7 yr). 

Recently, Weidenschilling (2003) has argued that metallicity enhancement is unlikely 
because Epstein-type drag, as described above, inadequately describes the forces on particles. 
Since the timescale for a particle to settle vertically toward the midplane is shorter than 
the timescale for a particle to drift radially by Epstein drag [by a factor of order rj ~ 
10~ 3 (r/AU) 1//2 ; see the next section], particles settle vertically into states of marginal Kelvin- 
Helmholtz turbulence (see Sekiya 1998) before drifting appreciably in the radial direction. 
Such turbulence exerts additional stresses on particles, stresses that are not accounted for by 
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Epstein drag; radial drift rates could, in principle, be altered significantly. Weidenschilling 
(2003) proceeds by assuming, as have other authors in different contexts (e.g., Goldreich & 
Ward 1973; Goodman & Pindor 2000, hereafter GP00), that this extra turbulent stress can 
be modeled as "plate drag": the turbulent stress on the particle layer is taken to be akin 
to the turbulent stress on a solid, rigid plate that is dragged through gas. The turbulent 
stresses so prescribed still induce inward radial drift of solids relative to gas, but the new, 
equilibrium drift velocities scale with stellocentric distance in such a way that particle pile- 
ups and metallicity enhancements do not occur in the MMSN. 

GP00 also employ the plate drag prescription to study the local, two dimensional dy- 
namics of the turbulent dust sheet. They discover an instability that causes in-plane per- 
turbations to grow over the local dynamical time. This instability was not recovered by 
Weidenschilling (2003) because the latter solved for the equilibrium state, in which unbal- 
anced gravitational forces on the dust layer match exactly the drag force; perturbations in 
radial momentum about this state were not considered. 1 

Our purpose here is to discard the usual plate drag prescription in favor of a self- 
consistent calculation of turbulent stresses based on the structure of the stratified particle 
layer. The essence of our approach is to solve for the turbulent diffusivity necessary to keep 
solids mixed vertically in the Kelvin-Helmholtz turbulent state, and then to employ this 
diffusivity in calculating rates of angular momentum transport. We find that the inclusion 
of turbulent stresses promotes local metal enrichment as particles accrete toward the central 
star; particles pile up radially more quickly with such stresses than without them. 

Our analysis below is restricted to millimeter-sized particles, though it can be extended 
to govern smaller sizes. As long as the Schmidt numbers of individual particles are nearly 
unity (see §2.2.2), the particles are well entrained in the gas and we may model the gas and 
particles as a single fluid. One millimeter is an interesting size regime from a number of 
perspectives. First, insofar as we would like to form planetesimals from seed particles that 
are as small as possible, with minimum recourse to collisional sticking, millimeter (or smaller) 
scales are more interesting to us than the meter scales that concern other works such as GP00. 
Second, collisional sticking might stall at sizes of a few millimeters, since relative velocities 
between particles tend to increase with increasing particle size, and relative velocities that are 
too high result in shattering of particle aggregates rather than sticking. Current consensus 
tells us that the threshold velocity that divides sticking from shattering is ~ 1 m s" 1 and 
that the corresponding maximum particle size generated by particle-particle collisions is ~ 1 
cm (Weidenschilling & Cuzzi 1993; Blum & Wurm 2000; Wurm, Blum, k Colwell 2001; 



1 These equilibrium states are referred to as "constant states" by GP00. 



-4- 



Chiang 2003). Third, modelling of millimeter- wave spectra of T Tauri and Herbig Ae star- 
disk systems offers evidence in favor of millimeter-sized grains dominating the particle mass 
near disk midplanes; the evidence is critically reviewed by Chiang (2003). Fourth, millimeter 
sizes characterize chondrules, the dominant constituent of the most primitive meteorites. 

In §2 we describe aerodynamic drift mechanisms for solids. Epstein drag is recapitulated 
in §2.1, while turbulent stresses are considered in §2.2. We review the plate drag prescription 
in §2.2.1 and introduce our new, self-consistent approach for calculating the transport of 
angular momentum by Kelvin-Helmholtz turbulence in §2.2.2. In §3 we apply our theory to 
a global simulation of the evolution of S p and demonstrate how particles pile up robustly. In 
§4 we argue more pointedly why our treatment of turbulent drag represents an improvement 
over the usual plate drag approximation. Concluding remarks are made in §5. Related 
discussions of the character and strength of turbulence in planetesimal forming disks are 
relegated to the appendices. 



2. Mechanisms for Radial Drift of Solids 



2.1. Epstein Drag 

Gas in protoplanetary disks rotates at speeds lower than the Keplerian speed, Vk, by 
an amount ^r/vx, where 

dPg/dlnr ^ ( c g 

2 Pg V K W 

is a dimensionless measure of radial pressure support and P g is the gas pressure. Typical 
models of the MMSN are characterized by i] ~ 10 _3 (r/AU) 1 / 2 . 

The Epstein aerodynamic drag law [equation (1)] implies that the stopping time of a 
particle of mass m p moving relative to gas is 




m p v rc i _ p s a 

£stop — n — • [o) 
^Ep PgCg 

For a particle of internal density p s = 3 g/cm 3 and radius a = 1 mm, the stopping time is 
likely shorter than the disk rotation period, r/vx = 1/^; to wit, 

^stop ~ KT 4 (^j)" , (4) 
where the index p describes the fall-off of gas surface density with distance, 
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stop *C 1, gas and solids rotate at nearly — to order fj(flt stop ) 
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velocity, 




(6) 



(Nakagawa et al. 1986). We make explicit the dependence of particle mass density, p p , on 
vertical height, z, above the midplane because vertical stratification of solids provides the 
only significant source of vertical shear, dv^/dz. This shear, in turn, drives Kelvin-Helmholtz 
turbulence. 

The radial drift speed of an individual particle can be derived via force balance in the 
radial and azimuthal directions (see, e.g., Goldreich & Ward 1973). Under the approximation 
that the azimuthal velocity of the particle equals the azimuthal velocity of the gas, the 
equation for radial force balance yields for the radial drift speed: 



The convention in this paper is that radially inward velocities are positive. Corrections to 
this expression are higher order in Qt stop . When the inertia of a collection of solids is taken 
into account, equation (7) becomes 



which is always smaller than t>Ep,md because p = p p + p g (Nakagawa et al. 1986). We refer 
to this dampening of the drift rate due to the collective inertia of the solids as the "inertial 
slow-down" effect. Note that while we have subscripted our velocities in equations (7) and 
(8) by "Ep" to denote the Epstein drag law, the form of the right-hand-side of equation (7) 
and the correction factor of pj p in equation (8) are independent of the drag law employed 
(e.g., for a > A, the appropriate drag law is due to Stokes: D oc a\c g v TC \). The same 
independence does not apply to the explicit form of tstop, which is given for the specific case 
of Epstein drag on the far right-hand side of equation (3). 

Particle pile-ups occur if the mass accretion rate in particles decreases with decreasing 
r at fixed time. Define 



^Ep,ind = ZV^tstopVK ■ 



(7) 




(8) 



S p oc r 



— n 



(9) 
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to be the surface density of solids. Further define the radial drift rate of solids to be, in 
general, 

v r ocr d . (10) 
Then the mass accretion rate in particles scales as 



Ti p rv r oc r , (11) 

where E — d — n + 1. If E > 0, then accretion results in particle pile-ups and local metal 
enrichment. 

For the specific case of Epstein drag, 

d Ep = ^-q + P, (12) 
3 

E Ep = -+p-n-q, (13) 

where the small inertial slow-down correction has been neglected. Here q describes the fall-off 
of midplane gas temperature with distance, 



T(xr~ q . (14) 

Reasonable values of q (see, e.g., Chiang & Goldreich 1997) span the range of 0.43-0.75. If 
the disk begins in a well-mixed state for which p = n, then E Ep > 1 initially and particles 
pile up as they accrete. Pile-ups cannot occur everywhere for all time. When a realistic disk 
with finite outer radius is considered, £ p at a fixed r will grow until particles at the outer 
radius drift past r. 

Figure 1 portrays the particle flux from Epstein drag, 

/Ep = PpVEp , (15) 

as a function of z for particles having a = 1 mm and p s = 3 g/cm 3 at r = 1 AU in Hayashi's 
(1981) model of the MMSN. We refer to this background nebular model hereafter as "model 
H." See YS02 for details of all nebular models used in this paper. Figure 2 is analogous to 
Figure 1, except that S p /E g is set to near the saturation limit, S p /S g = 0.99£ PiCr /£ g ; for 
the case of model H at r = 1 AU, this represents a factor of ~17 enhancement over solar 
metallicity. Distributions of particle density with height are computed using the model of 
Sekiya (1998) for which the Richardson number equals 1/4 everywhere. The fluxes are scaled 
to /0 g i>Ep,ind) which is independent of p p . 
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Fig. 1. — Radial particle fluxes from turbulent stresses, Epstein drag, and their sum vs. height 
above the midplane for millimeter-sized particles in model H at r ~ 1 AU. The vertically 
integrated abundance of chondrule-like particles is fixed at its cosmic value. Positive fluxes 
represent inflow. 



Before discussing the turbulent contribution to radial drift rates of solids, we prove that 
vertical settling of particles into the Kelvin-Helmholtz turbulent state occurs well before 
particles drift radially inward by Epstein drag. The timescale for vertical settling equals 

H 1 

taBtb ^Q2^f ' (16) 

X.O s top ^ L ''Stop 

where H g = c g /Q equals the vertical scale height of the gas, Vt 2 z equals the vertical compo- 
nent of stellar gravity, and we have set z = H g in the last equality in working to order-of- 
magnitude. The timescale for radial drift equals 

t Ep ~ — • (17) 
Substituting (7) and (8) into (17) reveals immediately that t se tt/^Ep ~ v(Pg/ P) 2 ^ 1- 
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2.2. Turbulent Stresses in Stratified Disks 

Forces exerted by gas turbulence are often modeled using the symmetric pressure tensor, 
P. In dynamical equilibrium, the advective transport of angular momentum balances the 
torque exerted by turbulent pressure fluctuations; to wit, 

dnr 2 r ( 1 dr 2 P r4> 1 8P H dP z<p 



Vr,tnih dr p \r 2 dr r <90 dz ) 

where i> r ,turb is the radial accretion velocity induced by turbulent stresses and (r, 0, z) are 
the usual cylindrical coordinates. Again, t> r ,turb > for radially inward accretion. The term 
proportional to P r( p dominates in active accretion disks that are rendered turbulent by, e.g., 
magneto-rotational instability (Balbus & Hawley 1991) or gravitational instability (Gammie 
2001). The disks considered here derive their turbulence from vertical shear instabilities, so 
that we set P r( f, = (see the appendices for quantitative justification). We further assume 
that the disk is axisymmetric so that dP^/dcp = 0. Equation (18) then reduces to 



^ r ,turb dVLr 2 ^ 1 dP z<j> 
r dr p dz 



(19) 



How can one model P Z( p in turbulent particle disks? In the next two subsections, we 
examine two prescriptions: one based on the conventional "plate drag" Ansatz, and a second 
rooted in the physics of stratified particle layers. 



2.2.1. The Plate Drag Approximation 

We review the "plate drag" treatment of turbulent stresses to connect to previous works 
[e.g., Goldreich & Ward (1973); GP00; Weidenschilling (2003)] and to contrast it with the 
new approach we champion in the next section §2.2.2. Section 4 criticizes the plate drag 
approximation more directly. 

Users of the plate drag prescription do not resolve the stress, P Z( p, as a function of z. 
Instead, P^ is vertically averaged and utilized to calculate an average accretion velocity. 
Averaging equation (19) over z yields 

for.turb) dttr 2 _ S 

r dr ~£(tf p )' 1 j 

where () denotes a density- weighted average, S(if p ) is the surface density of solids and gas 
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within the particle layer of half-thickness if p , and 

Hp ^fdz = 2P z 4z = H p ). (21) 

We have used the fact that P z< j, vanishes for \z\> H p and is an odd function of z. The heart 
of the plate drag prescription lies in setting the drag force per unit area, S, equal to 

where Av^ ~ tjvk is the characteristic azimuthal velocity difference between a dense particle- 
dominated midplane and overlying particle-free gas. 

Prescription (22) is useful only to the extent that Re* and any radial variations in this 
quantity are known. Goldreich & Ward (1973) assume that Re* = 500, citing an analysis 
of tidal boundary layers on the Earth's ocean floor. More recently, Cuzzi et al. (1993) and 
Dobrovolskis et al. (1999) obtain values of Re* ranging anywhere from 20 to 180 by using 
(the square of) three experimentally determined boundary lengths. Whether the results of 
these terrestrial-scale experiments can be extrapolated to astronomical settings has not been 
rigorously demonstrated. 

Insertion of (22) into (20) yields the equilibrium accretion velocity 

(Ur.turb) = Opiate = g R ^ , (23) 

where we have replaced H(H P ) with X p for consistency with other authors (although this 
replacement is strictly valid only if p p ^> p g ) and for ease when making analytic scalings. 
Equation (23) implies that the indices d and E introduced in §2.1 take the values <i p i a te = 
1 — 3g/2 — p + n and E p \ ate = 2 — 3q/2 — p, respectively, if we assume that Re* is spatially 
constant. 

In model H, p = n = 3/2 and q = 1/2. Then E p i atc = —1/4; solids drain from the 
disk onto the star and do not pile up. This is the result obtained by Weidenschilling (2003), 
who went on to show that adding the particle flux due to Epstein drag to that due to plate 
drag caused metal enhancement and depletion effects to cancel nearly exactly. We offer two 
comments regarding these results. The first is that whether the equilibrium solutions yield 
particle pile-ups is highly model-dependent. For example, an initial surface density profile 
for which p — n — 1 and q — 1/2 gives -E p i a te = 1/4. 2 Our second comment is that such 



2 For particles larger than the gas mean free path, the drag law is due to Stokes (see §1), and E = 



-q/2-n, 
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Fig. 2. — Same as Figure 1, except the surface density of solids has been increased to near 
the saturation limit: S p = 0.99E pcr . 

equilibrium solutions neglect the instability discovered by GP00, which yields metallicity 
enhancement on dynamical timescales when the plate drag prescription applies. 

In the next section we derive an alternative prescription for turbulent stresses based on 
the structure of Kelvin- Helmholtz turbulent flows. 



2.2.2. A Self- Consistent Description of Turbulent Stresses 
Proceeding from equation (19), we express the turbulent stress, 

i^ = P^, (24) 

in terms of a momentum diffusivity (viscosity), u z , which controls the efficiency of vertical 
mixing of azimuthal momentum. This same diffusivity also characterizes vertical mixing 
of particles; the ratio of the diffusivity of momentum to that of particles is the Schmidt 



i.e., no pile- up for decreasing temperature and density profiles. Loss of such large particles to accretion onto 
the star is another reason why we restrict ourselves to the millimeter sizes. 
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number, Sc; it is unity for the small particles considered here for which Qt stop <C 1 and 
that are entrained in all but the smallest turbulent eddies (Cuzzi et al. 1993). Given the 
equivalence of these diffusivities, we may derive v z from vertical profiles of particle density 
as computed by Sekiya (1998). Our procedure is developed as follows. 

The equilibrium distribution of particle density, p p (z), reflects a detailed balance be- 
tween turbulent diffusion and gravitational settling. The speed of upward particle diffusion 
is given by 

dhipp 

Wdiff = -v z dz , (25) 
while the downward vertical settling speed is nearly terminal, 

w sett = -Q 2 zt stop . (26) 

Vertical flux balance, tUdiff + ^sett = , allows us to calculate the turbulent viscosity in terms 
of known quantities: 

"-"Sts- < 27 > 

Note that the local viscosity increases with stopping time and with height above the mid- 
plane. More vigorous turbulence is required to keep larger particles aloft and at greater 
distances above the midplane where the vertical component of stellar gravity is stronger. In 
addition, strong gradients in particle density survive only in regions of low viscosity. 

We use (6) and (27) in (24) to calculate the relevant component of the stress tensor as 

Pz4> « -vrzn\ top ^p p . (28) 
Insertion of (28) into (19) yields the resulting equilibrium drift speed, 

V ^ h ~ -pn^T- ~ VEp ' ind 7^ U J • ( } 

Since the gas scale height, H g , is larger than the particle scale height (H p /H g ~ ^/rj 1), we 
have ignored vertical gradients in gas properties in writing equation (29). Note that f rj turb 
differs from the Epstein drag rate only via the shapes of the gas and particle density profiles. 
Since p p z/p is antisymmetric about the midplane and decreases near the top of the particle 
layer (as p p — > 0), equation (29) implies radial inflow near the midplane and radial outflow 
at greater heights. 

Note further that with the exception of the last step in equation (29), our treatment in 
this section does not depend on the validity of the Epstein drag law; the only requirement 
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is that the Schmidt number be unity, or equivalently that particles be small enough that 
f^stop <C 1- Strict adoption of the results in Cuzzi et al. (1993) — see, in particular, their 
Figure 2 — tells us that the Schmidt number for millimeter-sized particles in a Hayashi nebula 
varies from Sc = 1 inside r = 1 AU to Sc = 10 at r = 100 AU. For particle sizes a > 1 cm, the 
Schmidt number varies more strongly with r. Unfortunately, while it is true asymptotically 
that flt s top <C 1 implies Sc = 1, the precise values for Sc cited in Cuzzi et al. (1993) depend 
on poorly constrained dimensionless numbers; see their sections 2.2 and 2.3. We work, for 
simplicity, under the assumption of constant Sc ~ 1; despite the uncertainties involved, even 
this latter requirement could be relaxed by introducing a spatially variable Schmidt number. 

The equilibrium particle mass flux due to turbulent stresses, 

/turb = Pp^r,turb , (30) 

is plotted in Figures 1 and 2 for the cases of solar and super-solar metallicities, together 
with the previously discussed fluxes due to Epstein drag. In the case of solar metallicity, the 
turbulent flux increases the midplane flux by ~ 50%. The relative importance of turbulent 
fluxes increases with increasing metallicity; for the metal-rich case displayed in Figure 2, 
turbulent fluxes dominate near the midplane. 

It appears that turbulent stresses enhance the mass accretion rates due to Epstein drag. 
We proceed to execute global simulations of radial drift to investigate the possibility of 
particle pile-ups in the presence of both Epstein drag and turbulent stresses. 



3. Global Evolution 

We solve the mass continuity equation in the radial direction. Integrated over height, 
the continuity equation reads 

where the total particle flux equals 

F = F Ep + F tnrh = / (/ Ep + / turb ) dz . (32) 

As before, we assume an axially symmetric disk that is everywhere in a state of marginal 
Kelvin-Helmholtz turbulence. 

The integrals of flux over height can be computed analytically for the Sekiya models 
once the midplane density, p p (0), is known. One finds 



p g r]rv E p,md 



4 ( \x ~ 1 ) + ^ ' ~ 



l + ^o-3x+yX 2 ) /v (33) 
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Pg?7™Ep,ind 12 

where 



{ [Au 2 + x(Hx - 9)] Tl + 3 - 3 X ) - 2^ X 2 ] C) (34) 



+ (35) 
7 Pg /fi 2 , (36) 

(37) 

(38) 



Uq 




p g + p p (0) 


\ = 


i+^=i+ 


n = 


\Jx 2 ~u 2 o, 


c = 









Although cumbersome, these exact expressions permit greater accuracy and computational 
efficiency over direct numerical integration. 

Input parameters and boundary values for our first-order partial differential equation 
(31) are as follows. Particles have p s = 3 g/cm 3 and a = 1 mm. The gas disk obeys 
fixed power laws [equations (5) and (14)] with specific values for indices and normalizations 
taken from either model H or model Af (see YS02). The initial surface density profile for 
solids is scaled to that of gas but with an exponential cut-off at large radius, S p (r, t = 0) = 
£ g exp(-x 2 )/200, where x = r/(200 AU). The choice of £ p /£ g = 5 x 1(T 3 at x < 1 
assumes that only refractory material at solar abundances has condensed; inclusion of ices 
would enlarge £ p /£ g . Our outer boundary condition imposes a "lid" around the disk's 
outer edge, F(r out = 350 AU) = 0. We note that an inflow boundary condition, dF/dr = 
constant at r = r out , yields nearly identical results because of the exponential cutoff in S p . 
No boundary condition is required at the inner boundary, r in = 0.5 AU, where material 
leaves the grid. 



3.1. Numerical Techniques 

As in all advection problems, we must obey the Courant condition. The time step, At, 
must be smaller than the time for the flow, at speed v r , to pass from one grid point to the 
next; i.e., the Courant number Co = Atv r /Ar < 1 everywhere, where Ar is the grid spacing. 
The computational cost is dominated by regions with small Ar/v. Since v r oc r d roughly, the 
most efficient (non-adaptive) grid spacing is Ar oc r d so that we can fix Co = 0.5 everywhere. 
Such a grid, which is logarithmic only for d — 1, is generated by the recursion relation: 

r j+1 = r j {l + erf- 1 ), (39) 

where the index j labels the gridpoint and e is a constant based on the number of radial 
gridpoints, N r , and the desired range of r- values. This same grid works for all times because 
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we calculate Co using v r = i>E P ,ind, so that d depends only on time-independent gas properties; 
see equation (12). Evolutionary increases in particle density only act to decrease the flow 
speed, which poses no danger of violating the Courant condition. 

We found (see Figure 5) that setting N r = 1000 permitted fast integration with more 
than sufficient accuracy. The finite difference technique for the solution of (31) is explicit, 
upwind, forward in time, and first order (see Press et al. 1992, chap. 19.1). This simple 
technique proved convergent, stable, and accurate. 

Analytic evaluation of the vertically integrated fluxes [equations (33) and (34)] requires 
that we first obtain p p (z = 0) given the current value of E p at every radial gridpoint. To this 
end, we employ a non-linear root finder based on an optimized Miiller's method. The root- 
finding algorithm requires good initial guesses; it is the rate-limiting step in our simulation. 

We halt the simulation once £ p = £ PjCr anywhere. At this point, p p (0) becomes formally 
infinite, and gravitational instability ensues. One can imagine continuing the simulation 
by removing any excess surface density above S pcr ; the surface density of agglomerated 
planetesimals could then be tracked as a function of position and time. 

3.2. Results 

We now come to our main result, the evolution of S p (r, t) subject to Epstein drag and 
turbulent stresses when the vertical density profiles, p p (z), are supplied by models of Sekiya 
(1998). These results improve on the analytic models of YS02 who considered only Epstein 
drag and ignored both the inertial slow-down effect and turbulent stresses. 

Figure 3 showcases the evolution of particle surface density in model H under two 
cases: one in which only Epstein drag (corrected for inertial slow-down) is included, and 
a second in which both Epstein drag and turbulent transport of angular momentum are 
included. Pile-ups of surface density occur whether or not turbulent stresses are included in 
the calculation. Turbulent stresses accelerate the inflow and promote particle pile-ups. The 
simulation including turbulent stresses was halted after 5.3 x 10 5 yrs when S p > £ PjCr at 
r = 2.5 AU. To obtain saturation densities at larger stellocentric radii, one could start with 
a larger disk or a more metal-rich disk in which more solids could be provided by condensible 
ices. 

For model H, planetesimal formation is first triggered at the outer boundary or "cliff 
edge" of the shrinking disk of solids. This behavior is explained as follows. If the actual 
particle surface density flattens from its initial profile — i.e., if |dlnE p /dlnr| < \p\ — then 
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r (AU) 

Fig. 3. — Evolution of the surface density, E p , of millimeter-sized particles with time for 
model H using 10 3 gridpoints. Solid lines include both Epstein drag and turbulent stresses 
as sources of radial drift; dashed lines ignore turbulent stresses for comparison. In both 
cases, E p grows and its profile flattens with time as the particle disk contracts radially. 
Gravitational instability is first triggered when E p crosses E PjCr at the shrinking outer edge, 
here r ~ 2.5 AU. 

E p will first exceed E PjCr at the outer boundary of the accreting disk of particles. Whether 
the radial profile of particle surface density flattens depends on how the timescale for local 
density amplification, E p /E p , scales with r. Equation (31) reveals that E p /E p oc r/v r oc r l ~ d . 
If we ignore corrections to v r from the inertial slow-down effect and from turbulent stresses, 
then d = d,E P — p — 3/2 for model H. Hence, E p /E p oc r" 1 / 2 , the radial profile of E p flattens 
with time, and planetesimal formation is triggered first at the outer boundary of the particle 
disk. 

Figure 4 displays the evolution for nebular model Af, for which initial density profiles 
are shallower (p = n = 1) and temperature profiles are steeper (q = 0.63) than for model 
H. Both Epstein drag and turbulent stresses are included. Results that include only Epstein 
drag are similar and not shown in Figure 4 so as not to clutter it. This simulation was halted 
after 1.2 x 10 6 yrs when E p > E P)Cr at r = 0.5 AU, the inner boundary of the simulation 
domain. By contrast to model H, planetesimal formation begins at the inner edge of the 
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Fig. 4. — Similar to Figure 3, except for model Af (p — n — 1, q — 0.63). Epstein drag 
and turbulent stresses are included. The evolution proceeds relatively slowly for this model 
because the velocity index d ~ g?e p = 0.87 (v r oc r d ) is smaller than for model H (d^p = 1-5), 
so that drift velocities are smaller for this model than for model H outside r ~ 1 AU. Since 
<iEp < 1, S p (r) steepens and first exceeds S PjCr at the inner boundary of the simulation, 
r in = 0.5 AU. 

solid disk for model Af; the difference arises because d = 0.87 < 1 for the latter model, so 
that Sp/Sp oc r +0 - 13 and the radial profile of particle surface density steepens (see the above 
discussion for model H). 



To gain confidence in our numerical results and to aid in their interpretation, we ran two 
test cases that employ simplified prescriptions for the particle flux. Since these simple test 
cases do not depend on the detailed features of Sekiya's (1998) model, we can also gauge the 
degree to which our results are model-independent. These simple models allow S p > E pcr 
because they are constructed without regard to the physics of Kelvin-Helmholtz turbulence, 
and computationally they are less intensive. 



3.3. 



Test Cases 




Fig. 5. — Demonstration that our numerical simulation converges to the analytic result of 
YS02, even when sharp density peaks are present. Inertial slow-down effects and turbulent 
stresses are ignored in this test case; their restoration would smooth away the density peaks 
seen here. The disk is described by model H with chondrule-like particles. 

We first work in the limit of low particle density, p p <C p g , and set the particle flux 
F = F Epind = S p f E p jind . This expression was used by YS02 to obtain analytic solutions 
for the global evolution of particle surface density. We check our numerical code against 
these analytic solutions in Figure 5 at different grid resolutions. We chose model H because 
it gives rise to density cusps that provide more rigorous tests of numerical resolution than 
the smooth density profiles generated by model Af. The numerical solutions converge to 
the analytic result with increasing N r ; by N r = 10 4 , the analytic and numerical results are 
virtually indistinguishable. (In practice, for the more computationally intensive models that 
require root finding, we chose N r = 10 3 .) 

Our other test case assumes that p p is vertically constant within the particle layer. This 
"slab" approximation is a useful bridge between the low-density limit (F = F Ep ) and the 
full treatment based on the Sekiya profiles. The vertically integrated fluxes due to Epstein 
drag and turbulent stresses under the slab approximation read, respectively, 
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Fig. 6. — Similar to Figure 3, except p p is assumed to be vertically constant within the 
particle layer. This "slab approximation" captures the salient features of the full treatment 
based on Sekiya's (1998) models, including inertial slow-down, turbulent speed-up, and the 
smoothing away of density peaks seen in Figure 5. This simulation employed N r = 10 3 grid 
points. 



where a = £ p / (2p g H p ). For comparison, in the limit of low particle density, F Ep / (2p g if p t> Ep ) = 
a increases without bound with increasing E p . In the slab treatment, accounting for the effect 
of inertial slow-down saturates the flux to values below Fe p for a > 1. 

Figure 6 portrays the evolution of particle surface density in nebular model H under 
the slab approximation. The evolution is remarkably similar to that displayed in Figure 
3. In generating Figure 6, we have taken H p = y/Rir/r = r/r/2. In both Figures 3 and 6, 
the density peaks seen in Figure 5 are smoothed away in similar fashion. Smoothing occurs 
because of flux saturation in regions where a > 1. Under the slab approximation, turbulent 
stresses alter the evolution of particle density more markedly than under the full treatment; 
however, it remains the case qualitatively that turbulent stresses abet particle pile-ups. 

Note that the particle pile-up occurs relatively quickly in Figure 5, while the evolutions 
shown in both Figures 3 and 6 evince similar degrees of inertial slow-down. The slow-down is 




(41) 
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Fig. 7. — Vertically integrated, radial particle fluxes versus £ p based on Sekiya's (1998) 
profiles. Results are shown for gas drag with and without turbulent stresses for several 
values of Toomre's Q parameter for the gas, Q g . Instabilities can occur where <9F/<9£ p < 0. 

most pronounced under the full treatment using the Sekiya profiles. These differences arise 
from the fact that the integral in (32), when performed at fixed E p , diminishes either when 
H p decreases or when p p (z) becomes more inhomogeneous. 

One qualitative difference arises between using the slab approximation and the Sekiya 
models. In the former case, the radial particle flux monotonically increases with S p and 
asymptotes to a constant value. In the latter case, the flux attains a global maximum at 
S p pa 0.6E PiCr , if turbulent stresses and Epstein drag are included. The maximum occurs at 
lower E p if only Epstein drift is considered. This turnover of the flux is shown in Figure 7. At 
E p > 0.6E p cr , the flux decreases with increasing £ p . The magnitude of this decrease depends 
on the strength of self-gravity of the gas, as measured by the parameter Q g = A/(y/2nip). 
Smaller values of Q g indicate stronger self-gravity and cause a larger drop in the flux. The 
turn-over occurs because as E p approaches E p cr , H p decreases and engenders greater inertial 
slow-down. In addition, as E p increases, larger amounts of solids are swallowed into the cusp 
in density near the midplane, where the inertial slow-down effect is greatest. Both of these 
behaviors of p p (z) with E p can be seen in Figure 1 of YS02. 

The turn-over in F versus E p can spur instabilities. Imagine that <9F/<9E P < and 
that a region of high E p exists downstream of a region of low £ p . The low-density region 
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transmits a comparatively large flux to the high-density region, while the high-density region 
parts with a comparatively small flux. The high-density region therefore amplifies its density 
at the expense of regions yet further downstream. In this way, the instability propagates 
downstream. Indeed, we believe that this instability manifested itself in our simulations 
in the form of a "saw-tooth" where the density fluctuated from one radial grid point to 
the next. The saw-tooth typically appeared when S p /S PiCr > 0.8 — 0.9, in the regime where 
<9F/(9£p < 0. Attempts to resolve the instability at higher grid resolutions were unsuccessful; 
more study is needed. 3 While the shape of the F-E p relation that we have computed raises 
the possibility that the effective saturation criterion should instead read S p > 0.6S PjCr , such 
a refinement may not be important in practice. If aerodynamic drift (and other mechanisms) 
can raise S p to within a factor of 2 of S PjCr , they are likely to raise it all the way. 



4. On the Validity of Plate Drag 

In §2.2.2 and §3 we treated turbulent stresses without recourse to the plate drag ap- 
proximation. Since the plate drag formulae are commonly employed in the literature, we 
explain here why we believe that their use is problematic at best and inappropriate at worst. 

We have already mentioned in §2.2.1 the large range of values that the critical pa- 
rameter, Re*, has historically taken and the difficulty in determining its true form. More 
fundamentally, the plate drag approximation begs the question of whether we can think of 
a slurry of particles as a monolithic entity. GP00 propose that such collective behavior is 
possible if wakes around particles overlap, analogous to the practice of "drafting" in bicycle 
racing. In the present context of the Epstein drag regime, these wakes are disturbances 
having tiny length scales in the free molecular flow. If one assumes that each solid particle 
leaves a "footprint" on the surrounding fluid flow that extends a distance of order the gas 
mean free path, A, then one can show that for p p ~ p g the probability that another particle 
lies within this sphere of influence is ~10~ 8 . 

Despite the failure of free molecular wakes to overlap, collective effects among solids 
are not precluded. Indeed, the "inertial slow-down" correction to the particle drift velocity 
[equation (7)] tells us that particles are well aware of each other if their collective density, 
p p , approaches p g . This communication arises when one accounts for the back-reaction of 
particles dragging on gas. We find that this effect describes most of the collective behavior 
of solids, even when turbulent stresses, which are also collective in nature, are included. 



3 One could perform a local stability analysis analogous to that carried out by GP00. This would require 
restoring the time derivatives to our momentum equation (18). 
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Thus, plate drag overestimates turbulent stresses. Our final criticism concerns the as- 
sumption inherent in the plate drag approximation that the particle-dominated midplane 
behaves as a geometrically thin Ekman layer. We proceed to show that the turbulent dif- 
fusivity in such an Ekman layer is so large that the particle layer cannot be as thin as 
assumed. 

Ekman boundary layer theory applies to a medium having viscosity z/Ek and rotating 
at frequency Q, and assumes that other physics (e.g., buoyancy) is negligible. The height of 
the Ekman layer is 

#Ek « V^Ek/fi • (42) 

The viscosity, z/ Ek , which cannot be directly measured in protoplanetary disks, is defined in 
terms of i^Ek and the critical Reynolds number, 

Re* = H ^ VK , (43) 

^Ek 

where the characteristic velocity difference across the boundary layer is Av^ ~ tjvk- From 
equations (42) and (43), we derive the turbulent diffusivity of the layer, 

UEk ~ (Re*)2 ' (44) 

where all quantities, including the critical parameter Re*, are assumed to be known. The 
corresponding thickness of the Ekman layer is 

* Ek ~^. (45) 

Equations (42)-(45) are standard in the literature (e.g., Goldreich & Ward 1973; Dobrovolskis 
et al. 1999; GPOO). 

We compare H E k to the thickness of the particle layer that is stirred by turbulence char- 
acterized by z/Ek- Equation (27) with v z = u Ek yields, to order-of-magnitude, this thickness 4 : 

H Ek t . ^ / ^Ek ^ H Ek ^ 
Y ^ 2 ^stop -\J fitstop 



4 The assumption from § 2.2.2 that Sc w 1, true for r2t stop <C 1, is still in effect. 
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Whenever f2t st0 p <C 1, #Ek,stir is inconsistently larger than i^Ek- GPOO appreciate this point 
(see their equation [14]) and conclude that the plate drag formula yields self-consistent 
results for particles large enough (i.e., a ~ 10 cm at r ~ 1 AU) that Qt stop ~ 1. If, however, 
f^stop <S 1, then a layer as thin as will be stirred thicker until vertical shear becomes 
too weak to stir it further. This is the marginally Kelvin-Helmholtz turbulent state that we 
have employed throughout this paper. 

It should not be too surprising that the turbulent boundary layer as pictured within 
the plate drag approximation cannot, in general, be accurately described as an Ekman layer. 
The large Rossby number of the layer, 

Ro = -=^- « Re* > 1 , (47) 

"-nEk 

implies that inertia is more significant than rotation. Vertical shear and buoyancy are impor- 
tant ingredients missing from the traditional Ekman layer description. The Kelvin-Helmholtz 
turbulent layer that we have considered has Ro ~ 1. Recent work by Ishitsu & Sekiya (2003) 
suggests that centrifugal and Coriolis forces weaken vertical shear instabilities. This finding 
only helps to further pave the way for planetesimal formation by gravitational instability. 
See Appendix A for further discussion of rotational effects on particle drift rates. 

Finally, we comment that the Ekman layers in the Earth's atmosphere are not subject 
to the crisis described here. In the Earth's case, atmospheric shear is generated by thermal 
forcing, which is unaffected by boundary layer turbulence. By contrast, vertical shear in 
protoplanetary disks is generated by the inertia of particles within the boundary layer. 

In conclusion, in the small particle limit where Jli stop <C 1, the plate drag formula seems 
to apply only if particles were effectively glued onto a fixed plate! 



5. Summary 

We have investigated radial drift rates of solid particles within protoplanetary disks 
that derive their turbulence from vertical shear near their strongly stratified midplanes. 
Previous models by YS02 calculate drift rates according to mean flow differences between 
the sub-Keplerian gas and the more nearly Keplerian particles. We improve on their work 
by incorporating the effects of collective particle inertia [equation (7)] and the contribution 
to particle accretion from turbulent stresses [equation (29)]. These turbulent stresses are 
calculated in a physically self-consistent manner; turbulent particle diffusivities are derived 
from the detailed balance of vertically upward and downward particle fluxes within states of 
marginal Kelvin-Helmholtz turbulence [as computed by Sekiya (1998)]. Our results reinforce 
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the conclusions of YS02 that particle drifts lead to particle pile-ups. These pile-ups become 
sites of planetesimal formation by gravitational instability once the surface density of solids 
exceeds the saturation limit. We have explained why turbulent stresses cannot be modeled 
using the traditional plate drag approximation in the (most interesting) limit of small particle 
size; our alternative approach leads us to conclude that turbulent stresses hasten particle 
pile-ups. 

Inertial slow-down — the diminishing ability of gas to frictionally sap the angular momen- 
tum of particles as the particle mass density increases above the gas mass density — resulted 
in longer accretion timescales and a more even spreading of solids over larger regions of 
space, as compared to the results of YS02. Inclusion of this effect delays the particle pile-up, 
but only mildly. 

How the story of planetesimal formation unfolds in actual circumstellar disks will depend 
on a variety of factors in addition to aerodynamic drift of solids. Photoevaporation by 
stellar ultraviolet photons and stellar winds will strip gas from the disk surface and leave 
solids near the midplane. Furthermore, a distribution of particle sizes and geometries (fluffy 
vs. compact) will smooth particle surface density profiles. 

Is planetesimal formation triggered first near the outer or inner edges of accreting par- 
ticle disks? The answer bears directly on the architecture of planetary systems, on whether 
the mass in a planetary system is weighted toward small stellocentric distances (as they 
seem to be in extrasolar systems possessing Jovian-mass planets inside r ~ 1 AU) or larger 
stellocentric distances (as they are in our solar system). We have quantified the question 
of inside-out versus outside-in planet formation in terms of the index, d, which describes 
the inward radial drift speed of particles (v r oc r d ). If d > 1, then particle pile-ups occur 
first at the outer edge of the particle disk. For our power-law disks, d ~ d^ = p — q + 1/2, 
where p and q describe the variation of gas surface density and temperature with disk radius, 
respectively. Our knowledge of these quantities is informed by observations of circumstellar 
disks at mid-infrared to millimeter wavelengths. To date, g-values near 0.5 seem difficult to 
avoid (see, e.g., Chiang et al. 2001; D'Alessio et al. 2001), while estimates for p- values via 
dust continuum observations are traditionally hampered by lack of spatial resolution and 
ignorance concerning the dust opacity and the dust-to-gas ratio. Improvements in spatial 
resolution are promised by future detectors such as the Stratospheric Observatory for In- 
frared Astronomy (SOFIA) and the Atacama Large Millimeter Array (ALMA), while direct 
measurements of molecular hydrogen gas are possible via its rotational transitions at wave- 
lengths of 17 and 28 /xm (e.g., Thi et al. 1999) and its ro-vibrational lines at 2.1 /im (e.g., 
Bary et al. 2003). 

As noted by YS02, observations of our solar system's Kuiper Belt point with ever in- 
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creasing conviction to the existence of a hard outer edge to the classical Kuiper Belt at 
r :=» 48 AU (e.g., Allen et al. 2002, and references therein). Were it not for the coincidence 
of the edge with the location of the 2:1 mean-motion resonance with Neptune, we might 
simply ascribe this edge to a particle pile-up and the necessity of exceeding a critical thresh- 
old metallicity for planetesimal formation. We therefore offer a slightly more complicated 
scenario to explain these observations. A hard, primordial edge to the solar system formed 
at r < 48 AU as a consequence of the processes described in YS02 and this work. As 
Neptune was scattered (stochastically) outward by planetesimal encounters (Fernandez & 
Ip 1984; Hahn & Malhotra 1999), its exterior resonances captured bodies and carried them 
to greater heliocentric distances (Malhotra 1995; Chiang & Jordan 2002). In this way, the 
2:1 (strongest, outermost) Neptunian resonance combed the primordial edge outward to its 
present location. 

A. N. Y. acknowledges support from a National Science Foundation Graduate Research 
Fellowship. E.I.C. acknowledges support from Hubble Space Telescope Theory Grant HST- 
AR-09514.01-A and National Science Foundation Planetary Astronomy Grant AST-0205892. 
We thank the referee, Jeremy Goodman, for a thoughtful report that helped to improve the 
logic and presentation of this paper. 



A. Modifications to Turbulent Drift Rates by Rotation 

We have modeled Kelvin-Helmholtz turbulence as occurring in a vertical, Cartesian 
shear flow, without regard to the Coriolis force. This force deflects turbulent eddies in the 
z-<p plane into the r-0 plane and thereby introduces a non-zero P r( p stress into equation 
(18) for particle drift rates. Here we justify the neglect of this term, (l/r 2 )<9(r 2 P r(? <,)/<9r, as 
compared to dP z ^/dz. Since the Rossby number of the turbulent layer is of order unity, 

Ro = £*~£*~l, (Al) 

the magnitude of P r ^ is the same as that of P Z( p. An equivalent way of seeing this is to write 
down an expression for P r< ^ analogous to equation (24), 



Pr<p 

which is of the same order as 



(A2) 
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P# = pv z ^. pv z pv z il , (A3) 

az rjr 

since v r ~ Then it is clear that (1 / 'r 2 )d(r 2 P r <j>) j 'dr ~ P r ^jr is smaller than dP Z(j> /dz ~ 
P^/iJ p by ~# p /r ~ 77 < 1. 



B. Required Levels of Disk Passivity 

This appendix quantifies the extent to which circumstellar disks must be passive for the 
processes described in this paper to occur. Here we are considering "anomalous" sources 
of turbulence and angular momentum transport not arising from vertical shear in stratified 
particle layers (e.g., the magneto-rotational instability). We first ask under what conditions 
radial accretion velocities arising from an anomalous viscosity exceed the particle accretion 
velocities derived in this paper. For the latter we employ t>E P (since f r ,turb < ^Ep), and for 
the former we employ the usual a-prescription, 

v a ~ a r rjv K , (Bl) 

where a r parameterizes the strength of radial transport of angular momentum. Then one 
necessary criterion for disk passivity reads 

a r < ^nt stop . (B2) 

For chondrule-like particles at r ~ 1 AU in the MMSN, this criterion is satisfied for disks 
with a r 10~ 4 . It could be satisfied in more active disks for larger particles or at larger 
disk radii. 

A second criterion demands that the disk be sufficiently quiescent to permit vertical 
settling of solids. From equation (27), we estimate the diffusivity generated by Richardson 
turbulence as v z ~ (rjVK) 2 t s to P , where we have used H p ~ rjr. We express this diffusivity in 
terms of an a- viscosity as a z = v z VL/c 2 g , where the subscript z reminds us that turbulence 
need not be isotropic. Then this second criterion requires that any additional source of 
turbulence apart from vertical shear be characterized by 



« 2 < f^stop • (B3) 
For our power-law disks containing millimeter-sized particles, this criterion reads a z < 
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10 7 (r/AU) 1+p q . For model H, 1 + p — q = 2. The most stringent requirements on disk 
passivity lie at small stellocentric distances. 

We note that a z < 10~ 7 is still well satisfied by molecular viscosity, for which a r ~ a z ~ 
X/H g ~ 10- 12 at r ~ 1 AU. 

Does the anomalous viscosity ever shut off? Observations of T Tauri stars (and young 
brown dwarfs) indicate that infrared excesses indicating the presence of disks almost always 
correlate with Ha emission and other accretion signatures (Kenyon & Hartmann 1995; Liu 
et al. 2003). However, this oft-cited correlation has a considerable amount of scatter. Figure 
4 of Kenyon & Hartmann 1995 clearly shows that a sizable fraction of systems with infrared 
excesses have small Ha equivalent widths. These may represent disks that are quiescent 
enough to form planetesimals. We note that accretion signatures at optical-to-ultraviolet 
wavelengths pertain to activity in the immediate vicinity of the star; they may implicate 
accretion within ~30 stellar radii of the stellar surface, but they are not diagnostic of condi- 
tions at much greater stellocentric distances. Furthermore, evidence for significant vertical 
settling of dust in disk photospheres and the implied absence of vertical stirring in the up- 
permost scale heights is presented by Chiang et al. (2001), based on modeling mid-to-far 
infrared fluxes of Herbig Ae systems. 

In principle, it is possible for disks to exhibit a r ^> a z and to thereby reconcile modest 
levels of accretion with vertical settling of solids and particle pile-ups. For weakly viscous 
flows having small Rossby number, the Taylor- Proudman theorem constrains rotating flow 
to be two dimensional, parallel to the midplane, with no variation of velocity with height z 
(Pedlosky 1987, p. 43). Indeed, in both the Earth's atmosphere and the Earth's oceans, the 
ratios of horizontal to vertical viscosities are roughly estimated to be ~10 4 (Pedlosky 1987, 
p. 185). By analogy to our problem, we could simultaneously demand a z ~ 10 -8 so as to 
allow vertical settling and a r ~ 10 -4 so as to yield observable accretion rates of M ~ 10 -9 - 
lO _8 M /yr while not violating (B2) throughout most of the disk. Thus, mildly accreting 
class II protostellar systems could well be sites of planetesimal formation via gravitational 
instability. 
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